home *** CD-ROM | disk | FTP | other *** search
Text File | 1990-10-14 | 52.1 KB | 1,898 lines |
- Newsgroups: comp.sources.misc
- X-UNIX-From: craig@weedeater.math.yale.edu
- subject: v15i026: Graphics Gems, Part 4/5
- from: Craig Kolb <craig@weedeater.math.yale.edu>
- Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc)
-
- Posting-number: Volume 15, Issue 26
- Submitted-by: Craig Kolb <craig@weedeater.math.yale.edu>
- Archive-name: ggems/part04
-
- #! /bin/sh
- # This is a shell archive. Remove anything before this line, then unpack
- # it by saving it into a file and typing "sh file". To overwrite existing
- # files, type "sh file -c". You can also feed this as standard input via
- # unshar, or by typing "sh <file", e.g.. If this archive is complete, you
- # will see the following message at the end:
- # "End of archive 4 (of 5)."
- # Contents: 2DClip/cross.c AALines/AALines.c AALines/AATables.c
- # AAPolyScan.c ConcaveScan.c Forms.c Interleave.c Sturm/sturm.c
- # Wrapped by craig@weedeater on Fri Oct 12 15:53:14 1990
- PATH=/bin:/usr/bin:/usr/ucb ; export PATH
- if test -f '2DClip/cross.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'2DClip/cross.c'\"
- else
- echo shar: Extracting \"'2DClip/cross.c'\" \(6788 characters\)
- sed "s/^X//" >'2DClip/cross.c' <<'END_OF_FILE'
- X
- X/*
- X * file cross.c:
- X * calculate the intersections
- X */
- X#include <math.h>
- X#include "GraphicsGems.h"
- X#include "line.h"
- X#include "box.h"
- X/*
- X * cross_calc:
- X *
- X * PURPOSE
- X * calculate the intersections between the polygon
- X * stored in poly and the linesegment stored in l
- X * and put the intersections into psol.
- X *
- X * poly pointer to the structure containing the polygon
- X * l pointer to the structure containing the linesegment
- X * psol pointer to the pointer where intersections are stored
- X * nsol current number of intersections stored
- X * nsmax maximum storage in psol for intersections
- X * if nsol exceeds nsmax additional storage is allocated
- X *
- X */
- Xcross_calc(poly, l, psol, nsol, nsmax)
- XCONTOUR *poly;
- XSEGMENT *l;
- XCLIST **psol;
- Xshort *nsol, nsmax;
- X{
- X SEGMENT *p;
- X CLIST *sol;
- X double s;
- X long x, y, a, b, c;
- X int psort(), type;
- X
- X sol = *psol;
- X p = poly->_s;
- X do {
- X/*
- X * calculate the a, b and c coefficients and determine the
- X * type of intersection
- X */
- X
- X a = (p->_to._y - p->_from._y)*(l->_to._x - l->_from._x) -
- X (p->_to._x - p->_from._x)*(l->_to._y - l->_from._y);
- X b = (p->_from._x - l->_from._x)*(l->_to._y - l->_from._y) -
- X (p->_from._y - l->_from._y)*(l->_to._x - l->_from._x);
- X c = (p->_from._x - l->_from._x)*(p->_to._y - p->_from._y) -
- X (p->_from._y - l->_from._y)*(p->_to._x - p->_from._x);
- X if(a == 0)
- X type = (b == 0) ? COINCIDE : PARALLEL;
- X else {
- X if(a > 0) {
- X if((b >= 0 && b <= a) &&
- X (c >= 0 && c <= a))
- X type = CROSS;
- X else
- X type = NO_CROSS;
- X }
- X else {
- X if((b <= 0 && b >= a) &&
- X (c <= 0 && c >= a))
- X type = CROSS;
- X else
- X type = NO_CROSS;
- X }
- X }
- X/*
- X * process the interscetion found
- X */
- X switch(type) {
- X case NO_CROSS: case PARALLEL:
- X break;
- X
- X case CROSS:
- X if(b == a || c == a || c == 0)
- X break;
- X if(b == 0 &&
- X p_where(&(p->_prev->_from), &(p->_to), l) >= 0)
- X break;
- X s = (double)b/(double)a;
- X if(l->_from._x == l->_to._x)
- X x = l->_from._x;
- X else
- X x = p->_from._x +
- X (int)((p->_to._x - p->_from._x)*s);
- X if(l->_from._y == l->_to._y)
- X y = l->_from._y;
- X else
- X y = p->_from._y +
- X (int)((p->_to._y - p->_from._y)*s);
- X
- X if(*nsol == nsmax) {
- X nsmax *= 2;
- X *psol = sol = (CLIST *) realloc(sol, nsmax*sizeof(CLIST));
- X }
- X sol[*nsol]._p._x = x;
- X sol[*nsol]._p._y = y;
- X sol[*nsol]._type = STD;
- X *nsol += 1;
- X break;
- X
- X case COINCIDE:
- X if(p_where(&(p->_prev->_from),
- X &(p->_next->_to), l) > 0)
- X break;
- X if(l->_from._x != l->_to._x) {
- X if((max(l->_from._x, l->_to._x) <
- X min(p->_from._x, p->_to._x) ) ||
- X (min(l->_from._x, l->_to._x) >
- X max(p->_from._x, p->_to._x)) )
- X break;
- X if(min(l->_from._x, l->_to._x) <
- X min(p->_from._x, p->_to._x) ) {
- X if(*nsol == nsmax) {
- X nsmax *= 2;
- X *psol = sol = (CLIST *) realloc(sol, nsmax*sizeof(CLIST));
- X }
- X sol[*nsol]._p._x = p->_from._x;
- X sol[*nsol]._p._y = p->_from._y;
- X sol[*nsol]._type = DELAY;
- X *nsol += 1;
- X }
- X if(max(l->_from._x, l->_to._x) >
- X max(p->_from._x, p->_to._x) ) {
- X if(*nsol == nsmax) {
- X nsmax *= 2;
- X *psol = sol = (CLIST *) realloc(sol, nsmax*sizeof(CLIST));
- X }
- X sol[*nsol]._p._x = p->_to._x;
- X sol[*nsol]._p._y = p->_to._y;
- X sol[*nsol]._type = DELAY;
- X *nsol += 1;
- X }
- X }
- X else {
- X
- X if((max(l->_from._y, l->_to._y) <
- X min(p->_from._y, p->_to._y) ) ||
- X (min(l->_from._y, l->_to._y) >
- X max(p->_from._y, p->_to._y)) )
- X break;
- X if(min(l->_from._y, l->_to._y) <
- X min(p->_from._y, p->_to._y) ) {
- X if(*nsol == nsmax) {
- X nsmax *= 2;
- X *psol = sol = (CLIST *) realloc(sol, nsmax*sizeof(CLIST));
- X }
- X sol[*nsol]._p._x = p->_from._x;
- X sol[*nsol]._p._y = p->_from._y;
- X sol[*nsol]._type = DELAY;
- X *nsol += 1;
- X }
- X if(max(l->_from._y, l->_to._y) >
- X max(p->_from._y, p->_to._y) ) {
- X if(*nsol == nsmax) {
- X nsmax *= 2;
- X *psol = sol = (CLIST *) realloc(sol, nsmax*sizeof(CLIST));
- X }
- X sol[*nsol]._p._x = p->_to._x;
- X sol[*nsol]._p._y = p->_to._y;
- X sol[*nsol]._type = DELAY;
- X *nsol += 1;
- X }
- X }
- X break;
- X }
- X p = p->_next;
- X } while(p != poly->_s);
- X qsort(sol, *nsol, sizeof(CLIST), psort);
- X}
- X
- X
- X/*
- X * p_where
- X *
- X * PURPOSE
- X * determine position of point p1 and p2 relative to
- X * linesegment l.
- X * Return value
- X * < 0 p1 and p2 lie at different sides of l
- X * = 0 one of both points lie on l
- X * > 0 p1 and p2 lie at same side of l
- X *
- X * p1 pointer to coordinates of point
- X * p2 pointer to coordinates of point
- X * l pointer to linesegment
- X *
- X */
- Xp_where(p1, p2, l)
- XPOINT *p1, *p2;
- XSEGMENT *l;
- X{
- X long dx, dy, dx1, dx2, dy1, dy2, p_1, p_2;
- X
- X dx = l->_to._x - l->_from._x;
- X dy = l->_to._y - l->_from._y;
- X dx1 = p1->_x - l->_from._x;
- X dy1 = p1->_y - l->_from._y;
- X dx2 = p2->_x - l->_to._x;
- X dy2 = p2->_y - l->_to._y;
- X p_1 = (dx*dy1 - dy*dx1);
- X p_2 = (dx*dy2 - dy*dx2);
- X if(p_1 == 0 || p_2 == 0)
- X return(0);
- X else {
- X if((p_1 > 0 && p_2 < 0) || (p_1 < 0 && p_2 > 0))
- X return(-1);
- X else
- X return(1);
- X }
- X}
- X
- X
- X/*
- X * p_inside
- X *
- X * PURPOSE
- X * determine if the point stored in pt lies inside
- X * the polygon stored in p
- X * Return value:
- X * FALSE pt lies outside p
- X * TRUE pt lies inside p
- X *
- X * p pointer to the polygon
- X * pt pointer to the point
- X */
- Xboolean p_inside(p, pt)
- XCONTOUR *p;
- XPOINT *pt;
- X{
- X SEGMENT l, *sp;
- X CLIST *sol;
- X short nsol = 0, nsmax = 2, reduce = 0, i;
- X boolean on_contour(), odd;
- X
- X l._from._x = p->_minx-2;
- X l._from._y = pt->_y;
- X l._to._x = pt->_x;
- X l._to._y = pt->_y;
- X sol = (CLIST *) calloc(2, sizeof(CLIST));
- X cross_calc(p, &l, &sol, &nsol, nsmax);
- X for(i=0; i<nsol-1; i++)
- X if(sol[i]._type == DELAY && sol[i+1]._type == DELAY)
- X reduce++;
- X free(sol);
- X odd = (nsol - reduce) & 0x01;
- X return(odd ? !on_contour(p, pt) : FALSE);
- X}
- X
- X/*
- X * function used for sorting
- X */
- Xpsort(p1, p2)
- XCLIST *p1, *p2;
- X{
- X if(p1->_p._x != p2->_p._x)
- X return(p1->_p._x - p2->_p._x);
- X else
- X return(p1->_p._y - p2->_p._y);
- X}
- X
- X
- X/*
- X * on_contour
- X *
- X * PURPOSE
- X * determine if the point pt lies on the
- X * contour p.
- X * Return value
- X * TRUE point lies on contour
- X * FALSE point lies not on contour
- X *
- X * p pointer to the polygon structure
- X * pt pointer to the point
- X */
- Xboolean on_contour(p, pt)
- XCONTOUR *p;
- XPOINT *pt;
- X{
- X SEGMENT *sp;
- X long dx1, dy1, dx2, dy2;
- X
- X sp = p->_s;
- X do {
- X if((pt->_x >= min(sp->_from._x, sp->_to._x)) &&
- X (pt->_x <= max(sp->_from._x, sp->_to._x)) ) {
- X dx1 = pt->_x - sp->_from._x;
- X dx2 = sp->_to._x - pt->_x;
- X dy1 = pt->_y - sp->_from._y;
- X dy2 = sp->_to._y - pt->_y;
- X if(dy1*dx2 == dy2*dx1)
- X return(TRUE);
- X }
- X sp = sp->_next;
- X } while(sp != p->_s);
- X return(FALSE);
- X}
- X
- END_OF_FILE
- if test 6788 -ne `wc -c <'2DClip/cross.c'`; then
- echo shar: \"'2DClip/cross.c'\" unpacked with wrong size!
- fi
- # end of '2DClip/cross.c'
- fi
- if test -f 'AALines/AALines.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'AALines/AALines.c'\"
- else
- echo shar: Extracting \"'AALines/AALines.c'\" \(5139 characters\)
- sed "s/^X//" >'AALines/AALines.c' <<'END_OF_FILE'
- X/* FILENAME: AALines.c [revised 17 AUG 90]
- X
- X AUTHOR: Kelvin Thompson
- X
- X DESCRIPTION: Code to render an anti-aliased line, from
- X "Rendering Anti-Aliased Lines" in _Graphics_Gems_.
- X
- X This is almost exactly the code printed on pages 690-693
- X of _Graphics_Gems_. An overview of the code is on pages
- X 105-106.
- X
- X LINK WITH:
- X AALines.h -- Shared tables, symbols, etc.
- X AAMain.c -- Calling code for this subroutine.
- X AATables.c -- Initialize lookup tables.
- X*/
- X
- X
- X#include "AALines.h"
- X
- X#define SWAPVARS(v1,v2) ( temp=v1, v1=v2, v2=temp )
- X
- X#define FIXMUL(f1,f2) \
- X ( \
- X (((f1&0x0000ffff) * (f2&0x0000ffff)) >> 16) + \
- X ((f1&0xffff0000)>>16) * (f2&0x0000ffff) + \
- X ((f2&0xffff0000)>>16) * (f1&0x0000ffff) + \
- X ((f1&0xffff0000)>>16) * (f2&0xffff0000) \
- X )
- X
- X/* HARDWARE ASSUMPTIONS:
- X/* * 32-bit, signed ints
- X/* * 8-bit pixels, with initialized color table
- X/* * pixels are memory mapped in a rectangular fashion */
- X
- X/* FIXED-POINT DATA TYPE */
- X#ifndef FX_FRACBITS
- X typedef int FX;
- X# define FX_FRACBITS 16 /* bits of fraction in FX format */
- X# define FX_0 0 /* zero in fixed-point format */
- X#endif
- X
- X/* ASSUMED MACROS:
- X/* SWAPVARS(v1,v2) -- swaps the contents of two variables
- X/* PIXADDR(x,y) -- returns address of pixel at (x,y)
- X/* COVERAGE(FXdist) -- lookup macro for pixel coverage
- X given perpendicular distance; takes a fixed-point
- X integer and returns an integer in the range [0,255]
- X/* SQRTFUNC(FXval) -- lookup macro for sqrt(1/(1+FXval^2))
- X accepts and returns fixed-point numbers
- X/* FIXMUL(FX1,FX2) -- multiplies two fixed-point numbers
- X and returns the product as a fixed-point number */
- X
- X/* BLENDING FUNCTION:
- X/* 'cover' is coverage -- in the range [0,255]
- X/* 'back' is background color -- in the range [0,255] */
- X#define BLEND(cover,back) ((((255-(cover))*(back))>>8)+(cover))
- X
- X/* LINE DIRECTION bits and tables */
- X#define DIR_STEEP 1 /* set when abs(dy) > abs(dx) */
- X#define DIR_NEGY 2 /* set whey dy < 0 */
- X
- X
- X/* pixel increment values
- X/* -- assume PIXINC(dx,dy) is a macro such that:
- X/* PIXADDR(x0,y0) + PIXINC(dx,dy) = PIXADDR(x0+dx,y0+dy) */
- Xstatic int adj_pixinc[4] =
- X { PIXINC(1,0), PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1) };
- Xstatic int diag_pixinc[4] =
- X { PIXINC(1,1), PIXINC(1,1), PIXINC(1,-1), PIXINC(1,-1) };
- Xstatic int orth_pixinc[4] =
- X { PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1), PIXINC(1,0) };
- X
- X/* Global 'Pmax' is initialized elsewhere. It is the
- X "maximum perpendicular distance" -- the sum of half the
- X line width and the effective pixel radius -- in fixed format */
- XFX Pmax;
- X
- X
- X/*************** FUNCTION ANTI_LINE ***************/
- X
- Xvoid Anti_Line ( X1, Y1, X2, Y2 )
- Xint X1, Y1, X2, Y2;
- X{
- Xint Bvar, /* decision variable for Bresenham's */
- X Bainc, /* adjacent-increment for 'Bvar' */
- X Bdinc; /* diagonal-increment for 'Bvar' */
- XFX Pmid, /* perp distance at Bresenham's pixel */
- X Pnow, /* perp distance at current pixel (ortho loop) */
- X Painc, /* adjacent-increment for 'Pmid' */
- X Pdinc, /* diagonal-increment for 'Pmid' */
- X Poinc; /* orthogonal-increment for 'Pnow'--also equals 'k' */
- Xchar *mid_addr, /* pixel address for Bresenham's pixel */
- X *now_addr; /* pixel address for current pixel */
- Xint addr_ainc, /* adjacent pixel address offset */
- X addr_dinc, /* diagonal pixel address offset */
- X addr_oinc; /* orthogonal pixel address offset */
- Xint dx,dy,dir; /* direction and deltas */
- XFX slope; /* slope of line */
- Xint temp;
- X
- X/* rearrange ordering to force left-to-right */
- Xif ( X1 > X2 )
- X { SWAPVARS(X1,X2); SWAPVARS(Y1,Y2); }
- X
- X/* init deltas */
- Xdx = X2 - X1; /* guaranteed non-negative */
- Xdy = Y2 - Y1;
- X
- X
- X/* calculate direction (slope category) */
- Xdir = 0;
- Xif ( dy < 0 ) { dir |= DIR_NEGY; dy = -dy; }
- Xif ( dy > dx ) { dir |= DIR_STEEP; SWAPVARS(dx,dy); }
- X
- X/* init address stuff */
- Xmid_addr = PIXADDR(X1,Y1);
- Xaddr_ainc = adj_pixinc[dir];
- Xaddr_dinc = diag_pixinc[dir];
- Xaddr_oinc = orth_pixinc[dir];
- X
- X/* perpendicular measures */
- Xslope = (dy << FX_FRACBITS) / dx;
- XPoinc = SQRTFUNC( slope );
- XPainc = FIXMUL( slope, Poinc );
- XPdinc = Painc - Poinc;
- XPmid = FX_0;
- X
- X/* init Bresenham's */
- XBainc = dy << 1;
- XBdinc = (dy-dx) << 1;
- XBvar = Bainc - dx;
- X
- Xdo
- X {
- X /* do middle pixel */
- X *mid_addr = BLEND( COVERAGE(abs(Pmid)), *mid_addr );
- X
- X /* go up orthogonally */
- X for (
- X Pnow = Poinc-Pmid, now_addr = mid_addr+addr_oinc;
- X Pnow < Pmax;
- X Pnow += Poinc, now_addr += addr_oinc
- X )
- X *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
- X
- X /* go down orthogonally */
- X for (
- X Pnow = Poinc+Pmid, now_addr = mid_addr-addr_oinc;
- X Pnow < Pmax;
- X Pnow += Poinc, now_addr -= addr_oinc
- X )
- X *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
- X
- X
- X /* update Bresenham's */
- X if ( Bvar < 0 )
- X {
- X Bvar += Bainc;
- X mid_addr = (char *) ((int)mid_addr + addr_ainc);
- X Pmid += Painc;
- X }
- X else
- X {
- X Bvar += Bdinc;
- X mid_addr = (char *) ((int)mid_addr + addr_dinc);
- X Pmid += Pdinc;
- X }
- X
- X --dx;
- X } while ( dx >= 0 );
- X}
- END_OF_FILE
- if test 5139 -ne `wc -c <'AALines/AALines.c'`; then
- echo shar: \"'AALines/AALines.c'\" unpacked with wrong size!
- fi
- # end of 'AALines/AALines.c'
- fi
- if test -f 'AALines/AATables.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'AALines/AATables.c'\"
- else
- echo shar: Extracting \"'AALines/AATables.c'\" \(5772 characters\)
- sed "s/^X//" >'AALines/AATables.c' <<'END_OF_FILE'
- X/* FILENAME: AATables.c [revised 18 AUG 90]
- X
- X DESCRIPTION: Initialization of lookup tables and frame buffer
- X for anti-aliased line rendering demo.
- X
- X LINK WITH:
- X AALines.h -- Shared variables and symbols for renderer.
- X AAMain.c -- Calling routine for renderer.
- X AALines.c -- Anti-aliased line rendering code.
- X*/
- X
- X#include <math.h>
- X#include "AALines.h"
- X
- X/* programs in this file */
- Xextern void Anti_Init();
- Xstatic void Sqrt_Init();
- X
- X/* globals defined here */
- Xchar *fbuff;
- XUFX *sqrtfunc=0;
- Xint sqrtcells=1024;
- Xint sqrtshift;
- X
- X/* AA sizes */
- Xfloat line_r=1.0; /* line radius */
- Xfloat pix_r=SQRT_2; /* pixel radius */
- XFX *coverage=0;
- Xint covercells=128;
- Xint covershift;
- X
- X
- X
- X
- X/* ************************************
- X * *
- X * Anti_Init *
- X * *
- X ************************************
- X
- X DESCRIPTION: Initialize everything for the anti-aliased line
- X renderer in 'AALines.c': allocate the frame buffer,
- X set up lookup tables, etc.
- X
- X For hints about initializing the coverage table, see
- X "Area of Intersection: Circle and a Thick Line" on pages
- X 40-42 of _Graphics_Gems_.
- X*/
- X
- X#define FLOAT_TO_CELL(flt) ((int) ((flt) * 255.0 + 0.5))
- X#define MAXVAL_CELL 255
- X
- Xvoid Anti_Init ( )
- X{
- Xint *thiscell;
- Xdouble maxdist,nowdist,incdist;
- Xint tablebits,radbits;
- Xint tablecells;
- Xstatic int tablesize=0;
- Xdouble fnear,ffar,fcover;
- Xdouble half,invR,invpiRsq,invpi,Rsq;
- Xdouble sum_r;
- Xdouble inv_log_2;
- Xextern FX Pmax;
- X
- X/* alloc & init frame buffer */
- Xfbuff = (char *) malloc( xpix*ypix );
- X {
- X register int i;
- X for ( i=xpix*ypix-1; i>=0; --i ) fbuff[i] = 0;
- X }
- X
- X/* init */
- Xinv_log_2 = 1.0 / log( 2.0 );
- Xsum_r = line_r + pix_r;
- Xtablebits = (int) ( log((double)covercells) * inv_log_2 + 0.99 );
- Xradbits = (int) ( log((double)sum_r) * inv_log_2 ) + 1;
- Xcovershift = FX_FRACBITS - (tablebits-radbits);
- X
- X/* constants */
- Xhalf = 0.5;
- XinvR = 1.0 / pix_r;
- Xinvpi = 1.0 / PI;
- XinvpiRsq = invpi * invR * invR;
- XRsq = pix_r * pix_r;
- X#define FRACCOVER(d) (half - d*sqrt(Rsq-d*d)*invpiRsq - invpi*asin(d*invR))
- X
- X/* allocate table */
- XPmax = FLOAT_TO_FX(sum_r);
- XPmax >>= covershift;
- Xtablecells = Pmax + 2;
- XPmax <<= covershift;
- Xif ( coverage && tablecells > tablesize )
- X { free( coverage ); coverage = 0; tablesize = 0; }
- Xif ( coverage == 0 )
- X {
- X coverage = (FX *) malloc( tablecells * sizeof(int) );
- X tablesize = tablecells;
- X }
- X
- X/* init for fill loops */
- Xnowdist = 0.0;
- Xthiscell = coverage;
- Xincdist = sum_r / (double)(tablecells-2);
- X
- X/* fill fat portion */
- Xif ( pix_r <= line_r )
- X {
- X maxdist = line_r - pix_r;
- X for ( ;
- X nowdist <= maxdist;
- X nowdist += incdist, ++thiscell
- X )
- X {
- X *thiscell = MAXVAL_CELL;
- X }
- X }
- X
- X/* fill skinny portion */
- Xelse
- X {
- X /* loop till edge of line, or end of skinny, whichever comes first */
- X maxdist = pix_r - line_r;
- X if ( maxdist > line_r )
- X maxdist = line_r;
- X for ( ;
- X nowdist < maxdist;
- X nowdist += incdist, ++thiscell
- X )
- X {
- X fnear = line_r - nowdist;
- X ffar = line_r + nowdist;
- X fcover = 1.0 - FRACCOVER(fnear) - FRACCOVER(ffar);
- X *thiscell = FLOAT_TO_CELL(fcover);
- X }
- X
- X /* loop till end of skinny -- only run on super-skinny */
- X maxdist = pix_r - line_r;
- X for ( ;
- X nowdist < maxdist;
- X nowdist += incdist, ++thiscell
- X )
- X {
- X fnear = nowdist - line_r;
- X ffar = nowdist + line_r;
- X fcover = FRACCOVER(fnear) - FRACCOVER(ffar);
- X *thiscell = FLOAT_TO_CELL(fcover);
- X }
- X }
- X
- X/* loop till edge of line */
- Xmaxdist = line_r;
- Xfor ( ;
- X nowdist < maxdist;
- X nowdist += incdist, ++thiscell
- X )
- X {
- X fnear = line_r - nowdist;
- X fcover = 1.0 - FRACCOVER(fnear);
- X *thiscell = FLOAT_TO_CELL(fcover);
- X }
- X
- X/* loop till max separation */
- Xmaxdist = line_r + pix_r;
- Xfor ( ;
- X nowdist < maxdist;
- X nowdist += incdist, ++thiscell
- X )
- X {
- X fnear = nowdist - line_r;
- X fcover = FRACCOVER(fnear);
- X *thiscell = FLOAT_TO_CELL(fcover);
- X }
- X
- X/* finish off table */
- X*thiscell = FLOAT_TO_CELL(0.0);
- Xcoverage[tablecells-1] = FLOAT_TO_CELL(0.0);
- X
- XSqrt_Init();
- X}
- X
- X
- X
- X
- X/* *******************************
- X * *
- X * Sqrt_Init *
- X * *
- X *******************************
- X
- X DESCRIPTION: Initialize the lookup table for the function
- X sqrt(1/(1+x^2))). The table takes a shifted fixed-point
- X value as an index and returns a fixed-point value. Input
- X values are in the range [0,1] inclusive. The number of
- X cells in the table is a power of two plus one (the extra
- X cell provides an entry for an input value of exactly 1).
- X
- X GLOBALS:
- X sqrtcells -- Number of cells to use in the table (must
- X be set before calling this routine). This number is
- X rounded up to the nearest power of two (the global
- X variable itself is unchanged).
- X sqrtshift -- Bits to shift a fixed-point (FX) number
- X to generate a table index.
- X sqrtfunc -- Lookup table for the function.
- X*/
- X
- Xstatic void Sqrt_Init ( )
- X{
- XUFX *thiscell;
- Xdouble nowval,incval;
- Xint tablebits;
- Xint tablecells;
- Xdouble one;
- X
- X/* init */
- Xtablebits = (int) ( log((double)sqrtcells) / log(2.0) + 0.999 );
- Xtablecells = (1 << tablebits) + 1; /* one more that requested */
- Xsqrtshift = FX_FRACBITS - tablebits;
- Xone = 1.0;
- X
- X/* allocate table */
- Xif ( sqrtfunc )
- X free( sqrtfunc );
- Xsqrtfunc = (UFX *) malloc( tablecells * sizeof(int) );
- X
- X/* init for fill loop */
- Xincval = one / (double)(tablecells-1); /* a negative power of two */
- X
- Xfor (
- X nowval = 0.0, thiscell = sqrtfunc;
- X nowval < 1.0;
- X nowval += incval, ++thiscell
- X )
- X {
- X *thiscell = FLOAT_TO_FX( sqrt(one/(one+nowval*nowval)) );
- X }
- X
- Xsqrtfunc[tablecells-1] = FLOAT_TO_FX( sqrt(0.5) );
- X}
- END_OF_FILE
- if test 5772 -ne `wc -c <'AALines/AATables.c'`; then
- echo shar: \"'AALines/AATables.c'\" unpacked with wrong size!
- fi
- # end of 'AALines/AATables.c'
- fi
- if test -f 'AAPolyScan.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'AAPolyScan.c'\"
- else
- echo shar: Extracting \"'AAPolyScan.c'\" \(7519 characters\)
- sed "s/^X//" >'AAPolyScan.c' <<'END_OF_FILE'
- X/*
- XFast Anti-Aliasing Polygon Scan Converstion
- Xby Jack Morrison
- Xfrom "Graphics Gems", Academic Press, 1990
- X*/
- X
- X/*
- X* Anti-aliased polygon scan conversion by Jack Morrison
- X*
- X* This code renders a polygon, computing subpixel coverage at
- X* 8 times Y and 16 times X display resolution for anti-aliasing.
- X* One optimization left out for clarity is the use of incremental
- X* interpolations. X coordinate interpolation in particular can be
- X* with integers. See Dan Field's article in ACM Transactions on
- X* Graphics, January 1985 for a fast incremental interpolator.
- X*/
- X#include "GraphicsGems.h"
- X
- X#define SUBYRES 8 /* subpixel Y resolution per scanline */
- X#define SUBXRES 16 /* subpixel X resolution per pixel */
- X#define MAX_AREA (SUBYRES*SUBXRES)
- X#define MODRES(y) ((y) & 7) /*subpixel Y modulo */
- X#define MAX_X 0x7FFF /* subpixel X beyond right edge */
- X
- Xtypedef struct SurfaceStruct { /* object shading surface info */
- X int red, green, blue; /* color components */
- X } Surface;
- X/*
- X* In real life, SurfaceStruct will contain many more parameters as
- X* required by the shading and rendering programs, such as diffuse
- X* and specular factors, texture information, transparency, etc.
- X*/
- X
- Xtypedef struct VertexStruct { /* polygon vertex */
- X Vector3 model, world, /* geometric information */
- X normal, image;
- X int y; /* subpixel display coordinate */
- X } Vertex;
- X
- XVertex *Vleft, *VnextLeft; /* current left edge */
- XVertex *Vright, *VnextRight; /* current right edge */
- X
- Xstruct SubPixel { /* subpixel extents for scanline */
- X int xLeft, xRight;
- X } sp[SUBYRES];
- X
- Xint xLmin, xLmax; /* subpixel x extremes for scanline */
- Xint xRmax, xRmin; /* (for optimization shortcut) */
- X
- X/* Compute sub-pixel x coordinate for vertex */
- Xextern int screenX(/* Vertex *v */);
- X
- X/* Interpolate vertex information */
- Xextern void vLerp(/* double alpha, Vertex *Va, *Vb, *Vout */);
- X
- X/* Render polygon for one pixel, given coverage area */
- X/* and bitmask */
- Xextern void renderPixel(/* int x, y, Vertex *V,
- X int area, unsigned mask[],
- X Surface *object */);
- X
- X/*
- X * Render shaded polygon
- X */
- XdrawPolygon(polygon, numVertex, object)
- X Vertex polygon[]; /*clockwise clipped vertex list */
- X int numVertex; /*number of vertices in polygon */
- X
- X Surface *object; /* shading parms for this object */
- X{
- X Vertex *endPoly; /* end of polygon vertex list */
- X Vertex VscanLeft, VscanRight; /* interpolated vertices */ /* at scanline */
- X double aLeft, aRight; /* interpolation ratios */
- X struct SubPixel *sp_ptr; /* current subpixel info */
- X int xLeft, xNextLeft; /* subpixel coordinates for */
- X int xRight, xNextRight; /* active polygon edges */
- X int i,y;
- X
- X/* find vertex with minimum y (display coordinate) */
- XVleft = polygon;
- Xfor (i=1; i<numVertex; i++)
- X if (polygon[i].y < Vleft->y)
- X Vleft = &polygon[i];
- XendPoly = &polygon[numVertex-1];
- X
- X/* initialize scanning edges */
- XVright = VnextRight = VnextLeft = Vleft;
- X
- X/* prepare bottom of initial scanline - no coverage by polygon */
- Xfor (i=0; i<SUBYRES; i++)
- X sp[i].xLeft = sp[i].xRight = -1;
- XxLmin = xRmin = MAX_X;
- XxLmax = xRmax = -1;
- X
- X/* scan convert for each subpixel from bottom to top */
- Xfor (y+Vleft->y; ; y++) {
- X
- X while (y == VnextLeft->y) { /* reached next left vertex */
- X VnextLeft = (Vleft=VnextLeft) + 1; /* advance */
- X if (VnextLeft > endPoly) /* (wraparound) */
- X VnextLeft = polygon;
- X if (VnextLeft == Vright) /* all y's same? */
- X return; /* (null polygon) */
- X xLeft = screenX(Vleft);
- X xNextLeft = screenX(VnextLeft);
- X }
- X
- X while (y == VnextRight->y) { /*reached next right vertex */
- X VnextRight = (Vright=VnextRight) -1;
- X if (VnextRight < polygon) /* (wraparound) */
- X VnextRight = endPoly;
- X xRight = screenX(Vright);
- X xNextRight = screenX(VnextRight);
- X }
- X
- X if (y>VnextLeft->y || y>VnextRight->y) {
- X /* done, mark uncovered part of last scanline */
- X for (; MODRES(y); y++)
- X sp[MODRES(y)].xLeft = sp[MODRES(y)].xRight = -1;
- X renderScanline(Vleft, Vright, y/SUBYRES, object);
- X return;
- X }
- X
- X/*
- X * Interpolate sub-pixel x endpoints at this y,
- X * and update extremes for pixel coherence optimization
- X */
- X
- X sp_ptr = &sp[MODRES(y)];
- X aLeft = (double)(y - Vleft->y) / (VnextLeft->y - Vleft->y);
- X sp_ptr->xLeft = LERP(aLeft, xLeft, xNextLeft);
- X if (sp_ptr->xLeft < xLmin)
- X xLmin = sp_ptr->xLeft;
- X if (sp_ptr->xLeft > xLmax)
- X xLmax = sp_ptr->xLeft;
- X
- X aRight = (double)(y - Vright->y) / (VnextRight->y
- X - Vright->y);
- X sp_ptr->xRight = LERP(aRight, xRight, xNextRight);
- X if (sp_ptr->xRight < xRmin)
- X xRmin = sp_ptr->xRight;
- X if (sp_ptr->xRight > xRmax)
- X xRmax = sp_ptr->xRight;
- X
- X if (MODRES(y) == SUBYRES-1) { /* end of scanline */
- X /* interpolate edges to this scanline */
- X vLerp(aLeft, Vleft, VnextLeft, &VscanLeft);
- X vLerp(aRight, Vright, VnextRight, &VscanRight);
- X renderScanline(&VscanLeft, &VscanRight, y/SUBYRES, object);
- X xLmin = xRmin = MAX_X; /* reset extremes */
- X xLmax = xRmax = -1;
- X }
- X }
- X}
- X
- X/*
- X * Render one scanline of polygon
- X */
- X
- XrenderScanline(Vl, Vr, y, object)
- X Vertex *Vl, *Vr; /* polygon vertices interpolated */
- X /* at scanline */
- X int y; /* scanline coordinate */
- X Surface *object; /* shading parms for this object */
- X{
- X Vertex Vpixel; /*object info interpolated at one pixel */
- X unsigned mask[SUBYRES]; /*pixel coverage bitmask */
- X int x; /* leftmost subpixel of current pixel */
- X
- X for (x=SUBXRES*FLOOR(xLmin/SUBXRES); x<=xRmax; x+=SUBXRES) {
- X vLerp((double)(x-xLmin)/(xRmax-xLmin), Vl, Vr, &Vpixel);
- X computePixelMask(x, mask);
- X renderPixel(x/SUBXRES, y, &Vpixel,
- X /*computePixel*/Coverage(x), mask, object);
- X }
- X}
- X
- X/*
- X * Compute number of subpixels covered by polygon at current pixel
- X*/
- X/*computePixel*/Coverage(x)
- X int x; /* left subpixel of pixel */
- X{
- X int area; /* total covered area */
- X int partialArea; /* covered area for current subpixel y */
- X int xr = x+SUBXRES-1; /*right subpixel of pixel */
- X int y;
- X
- X /* shortcut for common case of fully covered pixel */
- X if (x>xLmax && x<xRmin)
- X return MAX_AREA;
- X
- X for (area=y=0; y<SUBYRES; y++) {
- X partialArea = MIN(sp[y].xRight, xr)
- X - MAX(sp[y].xLeft, x) + 1;
- X if (partialArea > 0)
- X area += partialArea;
- X }
- X return area;
- X}
- X
- X/* Compute bitmask indicating which subpixels are covered by
- X * polygon at current pixel. (Not all hidden-surface methods
- X * need this mask. )
- X*/
- XcomputePixelMask(x, mask)
- X int x; /* left subpixel of pixel */
- X unsigned mask[]; /* output bitmask */
- X{
- X static unsigned leftMaskTable[] =
- X { 0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF, 0x0FFF, 0x07FF, 0x03FF,
- X 0x01FF, 0x00FF, 0x007F, 0x003F, 0x001F, 0x000F, 0x0007,
- X 0x0003, 0x0001 };
- X static unsigned rightMaskTable[] =
- X { 0x8000, 0xC000, 0xE000, 0xF000, 0xF800, 0xFC00,
- X 0xFE00, 0xFF00, 0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
- X 0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF };
- X unsigned leftMask, rightMask; /* partial masks */
- X int xr = x+SUBXRES-1; /* right subpixel of pixel */
- X int y;
- X
- X/* shortcut for common case of fully covered pixel */
- X if (x>xLmax && x<xRmin) {
- X for (y=0; y<SUBYRES; y++)
- X mask[y] = 0xFFFF;
- X } else {
- X for (y=0; y<SUBYRES; y++) {
- X if (sp[y].xLeft < x) /* completely left of pixel*/
- X leftMask = 0xFFFF;
- X else if (sp[y].xLeft > xr) /* completely right */
- X leftMask = 0;
- X else
- X leftMask = leftMaskTable[sp[y].xLeft -x];
- X
- X if (sp[y].xRight > xr) /* completely */
- X /* right of pixel*/
- X rightMask = 0xFFFF;
- X else if (sp[y].xRight < x) /*completely left */
- X rightMask = 0;
- X else
- X rightMask = rightMaskTable[sp[y].xRight -x];
- X mask[y] = leftMask & rightMask;
- X }
- X }
- X}
- END_OF_FILE
- if test 7519 -ne `wc -c <'AAPolyScan.c'`; then
- echo shar: \"'AAPolyScan.c'\" unpacked with wrong size!
- fi
- # end of 'AAPolyScan.c'
- fi
- if test -f 'ConcaveScan.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'ConcaveScan.c'\"
- else
- echo shar: Extracting \"'ConcaveScan.c'\" \(4931 characters\)
- sed "s/^X//" >'ConcaveScan.c' <<'END_OF_FILE'
- X/*
- X * Concave Polygon Scan Conversion
- X * by Paul Heckbert
- X * from "Graphics Gems", Academic Press, 1990
- X */
- X
- X/*
- X * concave: scan convert nvert-sided concave non-simple polygon with vertices at
- X * (point[i].x, point[i].y) for i in [0..nvert-1] within the window win by
- X * calling spanproc for each visible span of pixels.
- X * Polygon can be clockwise or counterclockwise.
- X * Algorithm does uniform point sampling at pixel centers.
- X * Inside-outside test done by Jordan's rule: a point is considered inside if
- X * an emanating ray intersects the polygon an odd number of times.
- X * drawproc should fill in pixels from xl to xr inclusive on scanline y,
- X * e.g:
- X * drawproc(y, xl, xr)
- X * int y, xl, xr;
- X * {
- X * int x;
- X * for (x=xl; x<=xr; x++)
- X * pixel_write(x, y, pixelvalue);
- X * }
- X *
- X * Paul Heckbert 30 June 81, 18 Dec 89
- X */
- X
- X#include <stdio.h>
- X#include <math.h>
- X#include "GraphicsGems.h"
- X
- X#define ALLOC(ptr, type, n) ASSERT(ptr = (type *)malloc((n)*sizeof(type)))
- X
- Xtypedef struct { /* window: a discrete 2-D rectangle */
- X int x0, y0; /* xmin and ymin */
- X int x1, y1; /* xmax and ymax (inclusive) */
- X} Window;
- X
- Xtypedef struct { /* a polygon edge */
- X double x; /* x coordinate of edge's intersection with current scanline */
- X double dx; /* change in x with respect to y */
- X int i; /* edge number: edge i goes from pt[i] to pt[i+1] */
- X} Edge;
- X
- Xstatic int n; /* number of vertices */
- Xstatic Point2 *pt; /* vertices */
- X
- Xstatic int nact; /* number of active edges */
- Xstatic Edge *active; /* active edge list:edges crossing scanline y */
- X
- Xint compare_ind(), compare_active();
- X
- Xconcave(nvert, point, win, spanproc)
- Xint nvert; /* number of vertices */
- XPoint2 *point; /* vertices of polygon */
- XWindow *win; /* screen clipping window */
- Xvoid (*spanproc)(); /* called for each span of pixels */
- X{
- X int k, y0, y1, y, i, j, xl, xr;
- X int *ind; /* list of vertex indices, sorted by pt[ind[j]].y */
- X
- X n = nvert;
- X pt = point;
- X if (n<=0) return;
- X ALLOC(ind, int, n);
- X ALLOC(active, Edge, n);
- X
- X /* create y-sorted array of indices ind[k] into vertex list */
- X for (k=0; k<n; k++)
- X ind[k] = k;
- X qsort(ind, n, sizeof ind[0], compare_ind); /* sort ind by pt[ind[k]].y */
- X
- X nact = 0; /* start with empty active list */
- X k = 0; /* ind[k] is next vertex to process */
- X y0 = MAX(win->y0, ceil(pt[ind[0]].y-.5)); /* ymin of polygon */
- X y1 = MIN(win->y1, floor(pt[ind[n-1]].y-.5)); /* ymax of polygon */
- X
- X for (y=y0; y<=y1; y++) { /* step through scanlines */
- X /* scanline y is at y+.5 in continuous coordinates */
- X
- X /* check vertices between previous scanline and current one, if any */
- X for (; k<n && pt[ind[k]].y<=y+.5; k++) {
- X /* to simplify, if pt.y=y+.5, pretend it's above */
- X /* invariant: y-.5 < pt[i].y <= y+.5 */
- X i = ind[k];
- X /*
- X * insert or delete edges before and after vertex i (i-1 to i,
- X * and i to i+1) from active list if they cross scanline y
- X */
- X j = i>0 ? i-1 : n-1; /* vertex previous to i */
- X if (pt[j].y <= y-.5) /* old edge, remove from active list */
- X delete(j);
- X else if (pt[j].y > y+.5) /* new edge, add to active list */
- X insert(j, y);
- X j = i<n-1 ? i+1 : 0; /* vertex next after i */
- X if (pt[j].y <= y-.5) /* old edge, remove from active list */
- X delete(i);
- X else if (pt[j].y > y+.5) /* new edge, add to active list */
- X insert(i, y);
- X }
- X
- X /* sort active edge list by active[j].x */
- X qsort(active, nact, sizeof active[0], compare_active);
- X
- X /* draw horizontal segments for scanline y */
- X for (j=0; j<nact; j+=2) { /* draw horizontal segments */
- X /* span 'tween j & j+1 is inside, span tween j+1 & j+2 is outside */
- X xl = ceil(active[j].x-.5); /* left end of span */
- X if (xl<win->x0) xl = win->x0;
- X xr = floor(active[j+1].x-.5); /* right end of span */
- X if (xr>win->x1) xr = win->x1;
- X if (xl<=xr)
- X (*spanproc)(y, xl, xr); /* draw pixels in span */
- X active[j].x += active[j].dx; /* increment edge coords */
- X active[j+1].x += active[j+1].dx;
- X }
- X }
- X}
- X
- Xstatic delete(i) /* remove edge i from active list */
- Xint i;
- X{
- X int j;
- X
- X for (j=0; j<nact && active[j].i!=i; j++);
- X if (j>=nact) return; /* edge not in active list; happens at win->y0*/
- X nact--;
- X bcopy(&active[j+1], &active[j], (nact-j)*sizeof active[0]);
- X}
- X
- Xstatic insert(i, y) /* append edge i to end of active list */
- Xint i, y;
- X{
- X int j;
- X double dx;
- X Point2 *p, *q;
- X
- X j = i<n-1 ? i+1 : 0;
- X if (pt[i].y < pt[j].y) {p = &pt[i]; q = &pt[j];}
- X else {p = &pt[j]; q = &pt[i];}
- X /* initialize x position at intersection of edge with scanline y */
- X active[nact].dx = dx = (q->x-p->x)/(q->y-p->y);
- X active[nact].x = dx*(y+.5-p->y)+p->x;
- X active[nact].i = i;
- X nact++;
- X}
- X
- X/* comparison routines for qsort */
- Xcompare_ind(u, v) int *u, *v; {return pt[*u].y <= pt[*v].y ? -1 : 1;}
- Xcompare_active(u, v) Edge *u, *v; {return u->x <= v->x ? -1 : 1;}
- END_OF_FILE
- if test 4931 -ne `wc -c <'ConcaveScan.c'`; then
- echo shar: \"'ConcaveScan.c'\" unpacked with wrong size!
- fi
- # end of 'ConcaveScan.c'
- fi
- if test -f 'Forms.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'Forms.c'\"
- else
- echo shar: Extracting \"'Forms.c'\" \(5341 characters\)
- sed "s/^X//" >'Forms.c' <<'END_OF_FILE'
- X/*
- XForms, Vectors,and Transforms
- Xby Bob Wallis
- Xfrom "Graphics Gems", Academic Press, 1990
- X*/
- X
- X/*----------------------------------------------------------------------
- XThe main program below is set up to solve the Bezier subdivision problem
- Xin "Forms, Vectors, and Transforms". The subroutines are useful in
- Xsolving general problems which require manipulating matrices via exact
- Xinteger arithmetic. The intended application is validating or avoiding
- Xtedious algebraic calculations. As such, no thought was given to
- Xefficiency.
- X----------------------------------------------------------------------*/
- X#define ABS(x) ((x)>(0)? (x):(-x))
- X#define N 4 /* size of matrices to deal with */
- Xint M[N][N] = /* Bezier weights */
- X{
- X 1, 0, 0, 0,
- X -3, 3, 0, 0,
- X 3, -6, 3, 0,
- X -1, 3, -3, 1,
- X};
- Xint T[N][N] = /* re-parameterization xform for top half */
- X{
- X 1, -1, 1, -1,
- X 0, 2, -4, 6,
- X 0, 0, 4, -12,
- X 0, 0, 0, 8
- X};
- Xmain ()
- X{
- X int i,
- X j,
- X scale,
- X gcd,
- X C[N][N],
- X S[N][N],
- X Madj[N][N],
- X Tadj[N][N],
- X Mdet,
- X Tdet;
- X
- X
- X Tdet = adjoint (T, Tadj); /* inverse without division by */
- X Mdet = adjoint (M, Madj); /* determinant of T and M */
- X matmult (Madj, Tadj, C);
- X matmult (C, M, S); /* Madj*Tadj*M -> S */
- X scale = gcd = Mdet * Tdet; /* scale factors of both determinants */
- X for (i = 0; i < N; i++) /* find the greatest common */
- X { /* demoninator of S and determinants */
- X for (j = 0; j < N; j++)
- X gcd = Gcd (gcd, S[i][j]);
- X }
- X scale /= gcd; /* divide everything by gcd to get */
- X for (i = 0; i < N; i++) /* matrix and scale factor in lowest */
- X { /* integer terms possible */
- X for (j = 0; j < N; j++)
- X S[i][j] /= gcd;
- X }
- X printf ("scale factor = 1/%d ", scale);
- X print_mat ("M=", M, N); /* display the results */
- X print_mat ("T=", T, N);
- X print_mat ("S=", S, N); /* subdivision matrix */
- X exit (0);
- X}
- XGcd (a, b) /*returns greatest common demoninator */
- Xint a, /* of (a,b) */
- X b;
- X{
- X int i,
- X r;
- X a = ABS (a); /* force positive */
- X b = ABS (b);
- X if (a < b) /* exchange so that a >= b */
- X {
- X i = b;
- X b = a;
- X a = i;
- X }
- X if (b == 0)
- X return (a); /* finished */
- X r = a % b; /* remainder */
- X if (r == 0)
- X return (b); /* finished */
- X else /* recursive call */
- X return (Gcd (b, r));
- X}
- X
- Xadjoint (A, Aadj) /* returns determinant of A */
- Xint A[N][N], /* input matrix */
- X Aadj[N][N]; /* output = adjoint of A */
- X{ /* must have N >= 3 */
- X int i,
- X j,
- X I[N], /* arrays of row and column indices */
- X J[N],
- X Isub[N], /* sub-arrays of the above */
- X Jsub[N],
- X cofactor,
- X det;
- X if (N < 3)
- X {
- X printf ("must have N >= 3\n");
- X exit (1);
- X }
- X for (i = 0; i < N; i++)
- X { /* lookup tables to select a */
- X I[i] = i; /* particular subset of */
- X J[i] = i; /* rows and columns */
- X }
- X det = 0;
- X for (i = 0; i < N; i++)
- X { /* delete ith row */
- X subarray (I, Isub, N, i);
- X for (j = 0; j < N; j++)
- X { /* delete jth column */
- X subarray (J, Jsub, N, j);
- X cofactor = determinant (A, Isub, Jsub, N - 1,(i + j) & 1);
- X if (j == 0) /* use 0th column for det */
- X det += cofactor * A[i][0];
- X Aadj[j][i] = cofactor;
- X }
- X }
- X return (det);
- X}
- Xdeterminant (A, I, J, n, parity)/* actually gets a sub-determinant */
- Xint A[N][N], /* input = entire matrix */
- X I[N], /* row sub-array we want */
- X J[N], /* col sub-array we want */
- X parity, /* 1-> flip polarity */
- X n; /* # elements in subarrays */
- X
- X{
- X int i,
- X j,
- X det,
- X j_,
- X Jsub[N];
- X if (n <= 2) /* call ourselves till we get down to */
- X { /* a 2x2 matrix */
- X det =
- X (A[I[0]][J[0]] * A[I[1]][J[1]]) -
- X (A[I[1]][J[0]] * A[I[0]][J[1]]);
- X if (parity)
- X det = -det;
- X return (det);
- X } /* if (n <= 2) */
- X det = 0; /* n > 2; call recursively */
- X i = I[0]; /* strike out 0th row */
- X for (j_ = 0; j_ < n; j_++)
- X { /* strike out jth column */
- X subarray (J, Jsub, n, j_);
- X j = J[j_]; /* I + 1 => struck out 0th row */
- X det += A[i][j] * determinant (A, I + 1, Jsub, n - 1, j_ & 1);
- X }
- X if (parity)
- X det = -det;
- X return (det);
- X}
- Xsubarray (src, dest, n, k) /* strike out kth row/column */
- Xint *src, /* source array of n indices */
- X *dest, /* dest array formed by deleting k */
- X n,
- X k;
- X{
- X int i;
- X for (i = 0; i < n; i++, src++)
- X if (i != k) /* skip over k */
- X *dest++ = *src;
- X}
- X
- Xmatmult (A, B, C) /* C = A*B */
- Xint A[N][N],
- X B[N][N],
- X C[N][N];
- X{
- X int i,
- X j,
- X k,
- X sum;
- X for (i = 0; i < N; i++)
- X {
- X for (k = 0; k < N; k++)
- X {
- X sum = 0;
- X for (j = 0; j < N; j++)
- X sum += A[i][j] * B[j][k];
- X C[i][k] = sum;
- X }
- X }
- X}
- Xprint_mat (string, mat, n)
- Xchar *string;
- Xint mat[N][N],
- X n;
- X{
- X int i,
- X j;
- X printf ("%s\n", string);
- X for (i = 0; i < n; i++)
- X {
- X for (j = 0; j < n; j++)
- X printf (" %8ld", mat[i][j]);
- X printf ("\n");
- X }
- X}
- END_OF_FILE
- if test 5341 -ne `wc -c <'Forms.c'`; then
- echo shar: \"'Forms.c'\" unpacked with wrong size!
- fi
- # end of 'Forms.c'
- fi
- if test -f 'Interleave.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'Interleave.c'\"
- else
- echo shar: Extracting \"'Interleave.c'\" \(6497 characters\)
- sed "s/^X//" >'Interleave.c' <<'END_OF_FILE'
- X/*
- XBit Interleaving for Quad- or Octrees
- Xby Clifford A. Shaffer
- Xfrom "Graphics Gems", Academic Press, 1990
- X*/
- X
- X#include "GraphicsGems.h"
- X#define B_MAX_DEPTH 14 /* maximum depth allowed */
- X
- X/* byteval is the lookup table for coordinate interleaving. Given a
- X 4 bit portion of the (x, y) coordinates, return the bit interleaving.
- X Notice that this table looks like the order in which the pixels of
- X a 16 X 16 pixel image would be visited. */
- Xint byteval[16][16] =
- X { 0, 1, 4, 5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 85,
- X 2, 3, 6, 7, 18, 19, 22, 23, 66, 67, 70, 71, 82, 83, 86, 87,
- X 8, 9, 12, 13, 24, 25, 28, 29, 72, 73, 76, 77, 88, 89, 92, 93,
- X 10, 11, 14, 15, 26, 27, 30, 31, 74, 75, 78, 79, 90, 91, 94, 95,
- X 32, 33, 36, 37, 48, 49, 52, 53, 96, 97,100,101,112,113,116,117,
- X 34, 35, 38, 39, 50, 51, 54, 55, 98, 99,102,103,114,115,118,119,
- X 40, 41, 44, 45, 56, 57, 60, 61,104,105,108,109,120,121,124,125,
- X 42, 43, 46, 47, 58, 59, 62, 63,106,107,110,111,122,123,126,127,
- X 128,129,132,133,144,145,148,149,192,193,196,197,208,209,212,213,
- X 130,131,134,135,146,147,150,151,194,195,198,199,210,211,214,215,
- X 136,137,140,141,152,153,156,157,200,201,204,205,216,217,220,221,
- X 138,139,142,143,154,155,158,159,202,203,206,207,218,219,222,223,
- X 160,161,164,165,176,177,180,181,224,225,228,229,240,241,244,245,
- X 162,163,166,167,178,179,182,183,226,227,230,231,242,243,246,247,
- X 168,169,172,173,184,185,188,189,232,233,236,237,248,249,252,253,
- X 170,171,174,175,186,187,190,191,234,235,238,239,250,251,254,255};
- X
- X/* bytemask is the mask for byte interleaving - masks out the
- X non-significant bit positions. This is determined by the
- X depth of the node. For example, a node of depth 0 is at the root.
- X Thus, there are no branchs and no bits are significant.
- X The bottom 4 bits (the depth) are always retained.
- X Values are in octal notation. */
- Xint bytemask[B_MAX_DEPTH + 1] = {017,
- X 030000000017,036000000017,037400000017,037700000017,
- X 037760000017,037774000017,037777000017,037777600017,
- X 037777740017,037777770017,037777776017,037777777417,
- X 037777777717,037777777777};
- X
- X
- Xlong *interleave(addr, x, y, depth, max_depth)
- X/* Return the interleaved code for a quadtree node at depth depth
- Xwhose upper left hand corner has coordinates (x, y) in a tree with maximum
- Xdepth max_depth. This function receives and returns a
- Xpointer to addr, which is either a long interger or (more typically)
- Xan array of long integers whose first integer contains the
- Xinterleaved code. */
- X long *addr;
- X int max_depth, depth;
- X int x, y;
- X{
- X
- X/* Scale x, y values to be consistent with maximum coord size */
- X/* and depth of tree */
- X x <<= (B_MAX_DEPTH - max_depth);
- X y <<= (B_MAX_DEPTH - max_depth);
- X
- X/* calculate the bit interleaving of the x, y values that have now
- X been appropriately shifted, and place this interleave in the address
- X portion of addr. Note that the binary representations of x and y are
- X being processed from right to left */
- X
- X *addr = depth;
- X *addr |= byteval[y & 03][x & 03] << 4;
- X *addr |= byteval[(y >> 2) & 017][(x >> 2) & 017] << 8;
- X *addr |= byteval[(y >> 6) & 017][(x >> 6) & 017] << 16;
- X *addr |= byteval[(y >> 10) & 017][(x >> 10) & 017] << 24;
- X *addr &= bytemask[depth];
- X
- X/* if there were unused portions of the x and y addresses then */
- X/* the address was too large for the depth values given. */
- X/* Return address built */
- X return (addr);
- X}
- X
- X
- X
- X/* The next two arrays are used in calculating the (x, y) coordinates
- X of the upper left-hand corner of a node from its bit interleaved
- X address. Given an 8 bit number, the arrays return the effect of
- X removing every other bit (the y bits preceed the x bits). */
- X
- Xint xval[256] = { 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
- X 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
- X 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
- X 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
- X 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
- X 12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
- X 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
- X 12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
- X 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
- X 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
- X 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
- X 4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
- X 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
- X 12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
- X 8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
- X 12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15};
- X
- X
- Xint yval[256] = { 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
- X 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
- X 4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
- X 4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
- X 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
- X 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
- X 4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
- X 4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
- X 8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
- X 8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
- X 12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
- X 12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
- X 8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
- X 8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
- X 12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
- X 12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15};
- X
- X
- X
- X
- X
- Xint getx(addr, max_depth)
- X/* Return the x coordinate of the upper left hand corner of addr for a
- X tree with maximum depth max_depth. */
- X long *addr;
- X int max_depth;
- X{
- X register x;
- X
- X x = xval[(*addr >> 4) & 017];
- X x |= xval[(*addr >> 8) & 0377] << 2;
- X x |= xval[(*addr >> 16) & 0377] << 6;
- X x |= xval[(*addr >> 24) & 0377] << 10;
- X x >>= B_MAX_DEPTH - max_depth;
- X return (x);
- X}
- X
- X
- X
- Xint QKy(addr, max_depth)
- X/* Return the y coordinate of the upper left hand corner of addr for a
- X tree with maximum depth max_depth. */
- X
- X long *addr;
- X int max_depth;
- X{
- X register y;
- X
- X y = yval[(*addr >> 4) & 017];
- X y |= yval[(*addr >> 8) & 0377] << 2;
- X y |= yval[(*addr >> 16) & 0377] << 6;
- X y |= yval[(*addr >> 24) & 0377] << 10;
- X y >>= B_MAX_DEPTH - max_depth;
- X return (y);
- X}
- X
- Xint getdepth(addr)
- X/* Return the depth of the node. Simply return the bottom 4 bits. */
- X
- X long *addr;
- X{
- X
- X return(*addr & 017);
- X}
- END_OF_FILE
- if test 6497 -ne `wc -c <'Interleave.c'`; then
- echo shar: \"'Interleave.c'\" unpacked with wrong size!
- fi
- # end of 'Interleave.c'
- fi
- if test -f 'Sturm/sturm.c' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'Sturm/sturm.c'\"
- else
- echo shar: Extracting \"'Sturm/sturm.c'\" \(5198 characters\)
- sed "s/^X//" >'Sturm/sturm.c' <<'END_OF_FILE'
- X
- X/*
- X * sturm.c
- X *
- X * the functions to build and evaluate the Sturm sequence
- X */
- X#include <math.h>
- X#include <stdio.h>
- X#include "solve.h"
- X
- X/*
- X * modp
- X *
- X * calculates the modulus of u(x) / v(x) leaving it in r, it
- X * returns 0 if r(x) is a constant.
- X * note: this function assumes the leading coefficient of v
- X * is 1 or -1
- X */
- Xstatic int
- Xmodp(u, v, r)
- X poly *u, *v, *r;
- X{
- X int k, j;
- X double *nr, *end, *uc;
- X
- X nr = r->coef;
- X end = &u->coef[u->ord];
- X
- X uc = u->coef;
- X while (uc <= end)
- X *nr++ = *uc++;
- X
- X if (v->coef[v->ord] < 0.0) {
- X
- X
- X for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
- X r->coef[k] = -r->coef[k];
- X
- X for (k = u->ord - v->ord; k >= 0; k--)
- X for (j = v->ord + k - 1; j >= k; j--)
- X r->coef[j] = -r->coef[j] - r->coef[v->ord + k]
- X * v->coef[j - k];
- X } else {
- X for (k = u->ord - v->ord; k >= 0; k--)
- X for (j = v->ord + k - 1; j >= k; j--)
- X r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
- X }
- X
- X k = v->ord - 1;
- X while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
- X r->coef[k] = 0.0;
- X k--;
- X }
- X
- X r->ord = (k < 0) ? 0 : k;
- X
- X return(r->ord);
- X}
- X
- X/*
- X * buildsturm
- X *
- X * build up a sturm sequence for a polynomial in smat, returning
- X * the number of polynomials in the sequence
- X */
- Xint
- Xbuildsturm(ord, sseq)
- X int ord;
- X poly *sseq;
- X{
- X int i;
- X double f, *fp, *fc;
- X poly *sp;
- X
- X sseq[0].ord = ord;
- X sseq[1].ord = ord - 1;
- X
- X
- X /*
- X * calculate the derivative and normalise the leading
- X * coefficient.
- X */
- X f = fabs(sseq[0].coef[ord] * ord);
- X fp = sseq[1].coef;
- X fc = sseq[0].coef + 1;
- X for (i = 1; i <= ord; i++)
- X *fp++ = *fc++ * i / f;
- X
- X /*
- X * construct the rest of the Sturm sequence
- X */
- X for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
- X
- X /*
- X * reverse the sign and normalise
- X */
- X f = -fabs(sp->coef[sp->ord]);
- X for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
- X *fp /= f;
- X }
- X
- X sp->coef[0] = -sp->coef[0]; /* reverse the sign */
- X
- X return(sp - sseq);
- X}
- X
- X/*
- X * numroots
- X *
- X * return the number of distinct real roots of the polynomial
- X * described in sseq.
- X */
- Xint
- Xnumroots(np, sseq, atneg, atpos)
- X int np;
- X poly *sseq;
- X int *atneg, *atpos;
- X{
- X int atposinf, atneginf;
- X poly *s;
- X double f, lf;
- X
- X atposinf = atneginf = 0;
- X
- X
- X /*
- X * changes at positve infinity
- X */
- X lf = sseq[0].coef[sseq[0].ord];
- X
- X for (s = sseq + 1; s <= sseq + np; s++) {
- X f = s->coef[s->ord];
- X if (lf == 0.0 || lf * f < 0)
- X atposinf++;
- X lf = f;
- X }
- X
- X /*
- X * changes at negative infinity
- X */
- X if (sseq[0].ord & 1)
- X lf = -sseq[0].coef[sseq[0].ord];
- X else
- X lf = sseq[0].coef[sseq[0].ord];
- X
- X for (s = sseq + 1; s <= sseq + np; s++) {
- X if (s->ord & 1)
- X f = -s->coef[s->ord];
- X else
- X f = s->coef[s->ord];
- X if (lf == 0.0 || lf * f < 0)
- X atneginf++;
- X lf = f;
- X }
- X
- X *atneg = atneginf;
- X *atpos = atposinf;
- X
- X return(atneginf - atposinf);
- X}
- X
- X/*
- X * numchanges
- X *
- X * return the number of sign changes in the Sturm sequence in
- X * sseq at the value a.
- X */
- Xint
- Xnumchanges(np, sseq, a)
- X int np;
- X poly *sseq;
- X double a;
- X
- X{
- X int changes;
- X double f, lf;
- X poly *s;
- X
- X changes = 0;
- X
- X lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
- X
- X for (s = sseq + 1; s <= sseq + np; s++) {
- X f = evalpoly(s->ord, s->coef, a);
- X if (lf == 0.0 || lf * f < 0)
- X changes++;
- X lf = f;
- X }
- X
- X return(changes);
- X}
- X
- X/*
- X * sbisect
- X *
- X * uses a bisection based on the sturm sequence for the polynomial
- X * described in sseq to isolate intervals in which roots occur,
- X * the roots are returned in the roots array in order of magnitude.
- X */
- Xsbisect(np, sseq, min, max, atmin, atmax, roots)
- X int np;
- X poly *sseq;
- X double min, max;
- X int atmin, atmax;
- X double *roots;
- X{
- X double mid;
- X int n1, n2, its, atmid, nroot;
- X
- X if ((nroot = atmin - atmax) == 1) {
- X
- X /*
- X * first try a less expensive technique.
- X */
- X if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
- X return;
- X
- X
- X /*
- X * if we get here we have to evaluate the root the hard
- X * way by using the Sturm sequence.
- X */
- X for (its = 0; its < MAXIT; its++) {
- X mid = (min + max) / 2;
- X
- X atmid = numchanges(np, sseq, mid);
- X
- X if (fabs(mid) > RELERROR) {
- X if (fabs((max - min) / mid) < RELERROR) {
- X roots[0] = mid;
- X return;
- X }
- X } else if (fabs(max - min) < RELERROR) {
- X roots[0] = mid;
- X return;
- X }
- X
- X if ((atmin - atmid) == 0)
- X min = mid;
- X else
- X max = mid;
- X }
- X
- X if (its == MAXIT) {
- X fprintf(stderr, "sbisect: overflow min %f max %f\
- X diff %e nroot %d n1 %d n2 %d\n",
- X min, max, max - min, nroot, n1, n2);
- X roots[0] = mid;
- X }
- X
- X return;
- X }
- X
- X /*
- X * more than one root in the interval, we have to bisect...
- X */
- X for (its = 0; its < MAXIT; its++) {
- X
- X mid = (min + max) / 2;
- X
- X atmid = numchanges(np, sseq, mid);
- X
- X n1 = atmin - atmid;
- X n2 = atmid - atmax;
- X
- X
- X if (n1 != 0 && n2 != 0) {
- X sbisect(np, sseq, min, mid, atmin, atmid, roots);
- X sbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
- X break;
- X }
- X
- X if (n1 == 0)
- X min = mid;
- X else
- X max = mid;
- X }
- X
- X if (its == MAXIT) {
- X fprintf(stderr, "sbisect: roots too close together\n");
- X fprintf(stderr, "sbisect: overflow min %f max %f diff %e\
- X nroot %d n1 %d n2 %d\n",
- X min, max, max - min, nroot, n1, n2);
- X for (n1 = atmax; n1 < atmin; n1++)
- X roots[n1 - atmax] = mid;
- X }
- X}
- END_OF_FILE
- if test 5198 -ne `wc -c <'Sturm/sturm.c'`; then
- echo shar: \"'Sturm/sturm.c'\" unpacked with wrong size!
- fi
- # end of 'Sturm/sturm.c'
- fi
- echo shar: End of archive 4 \(of 5\).
- cp /dev/null ark4isdone
- MISSING=""
- for I in 1 2 3 4 5 ; do
- if test ! -f ark${I}isdone ; then
- MISSING="${MISSING} ${I}"
- fi
- done
- if test "${MISSING}" = "" ; then
- echo You have unpacked all 5 archives.
- rm -f ark[1-9]isdone
- else
- echo You still need to unpack the following archives:
- echo " " ${MISSING}
- fi
- ## End of shell archive.
- exit 0
-
-